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SUMMARY 


- — 'j  In  thisyffnal  Technical  Report  an  analagous  Bayesian  approach 
is  given  to  the  classical  methodology  for  calculating  lower  bounds  on 
system  reliability  as  formulated  in  the  Maximus  Handbook  of  February 


1980. ^ 

y 

^  For  ease  of  comparison  the  Maximus  format  is  adhered  to  as  far 
as  possible.  All  theoretical  work  is  given  in  appendices,  as  is  the 
description  of  the  computer  programs  that  were  developed  to  facilitate 
calculations.  ' 


^In  addition  to  the  pass/fail  test  data  considered  in  Maximus  the 
case  of  exponential  times  to  failure  of  components  is  treated. 


)/ 


vy >f*  so  . 

•  Key  Words:-  Bayek,  lower  limits,  system  reliability, 
component  test  data,  reduction  methods. 
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1.  INTRODUCTION 

It  is  pertinent  to  refer  to  remarks  made  in  the  Introduction  to 
the  Maximus  Handbook.  The  first  remark  recognises  the  need  for 
continued  research  and  analysis  and  that  the  methodology  of  Maximus 
should  be  considered  as  interim  in  nature.  The  second  remark  states 
that  the  limitation  is  to  Classical  Statistical  Theory  and  that 
extension  of  the  scope  of  the  work  to  cover  Bayesian  methods  is 
deferred . 

The  present  report  addresses  these  two  points.  Of  principal 
interest  is  the  Bayesian  analog  of  the  Classical  Maximus  methodology. 
However  during  the  course  of  development,  it  became  necessary  to  take 
a  fresh  look  at  the  latter  and  one  consequence  is  that  we  suggest 
improvements  to  the  Lindstrom  and  Madden  method  for  series  systems 
and  also  for  the  case  of  repeated  components. 

It  is  not  surprising  that  there  should  be  a  link  between 
Classical  and  Bayesian  approaches  to  reducing  component  information 
to  'equivalent'  system  information.  This  is  due  to  the  well  known 
relationship  between  partial  Binomial  sums  and  the  Incomplete  Beta 
function  (Chapter  2). 

The  Bayesian  Statistician  expresses  uncertainty  about  the 
unknown  component  reliabilities  through  prior  probability 
distributions.  For  any  given  component  the  test  data  is  used  to 
construct  a  likelihood  function  which  is  then  combined  with  the 
prior  probability  distribution,  using  Bayes'  theorem,  to  give  a 
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posterior  distribution  of  component  reliability.  The  system 
reliability  is  then  the  appropriate  function  (the  system  reliability 
function)  of  component  reliabilities.  The  determination  of  the 
exact  system  reliability  distribution  is  usually  formidably  difficult, 
hence  the  search  for  simple  but  good  approximate  methods,  including 
the  reduction  methods  developed  herein. 

The  intepretation  of  interval  estimate  is  different  using 
Classical  and  Bayes  approaches.  For  example,  lower  90%  limits 
have  the  following  interpretations. 

(i)  Classical :  There  is  at  least  a  0.90  probability  that  the 
lower  confidence  limit  (a  random  variable)  will  take  on  a 
value  lower  than  the  fixed  but  unknown  system  reliability. 

(ii)  Bayesian:  There  is  at  least  a  0.90  probability  that  the 

unknown  system  reliability  (a  random  variable)  is  greater 
than  the  lower  limit. 

When  random  variables  and  parameters  are  continuous  the  phrase 
'at  least'  can  be  omitted  in  (i)  and  also  in  (ii)  for  continuously 
distributed  parameters. 

Classical  limits  are  referred  to  as  confidence  limits.  Some¬ 
times  Bayesian  limits  are  referred  to  as  Bayesian  confidence  limits 
but  we  shall  reserve  the  use  of  the  word  confidence  for  the  classical 
case  and  simply  use  the  term  Bayesian  limits.  These  are  appropriate 
percentage  points  of  the  posterior  distribution  of  reliability. 

The  terms  credible  intervals  and  credible  limits  are  increasingly  used 
in  the  Bayesian  context. 
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Like  Maximus  we  restrict  our  attention  to  series /pa  rail  el  systems, 
for  other  structures  such  as  k  out  of  n  or  quorum  structures  the 
methods  advocated  in  DAJA37-82-C-0736  can  be  used.  These  methods, 
based  on  asymptotic  expansions,  can  also  be  used  to  provide  a  useful 
check  on  the  accuracy  of  the  reduction  methods  herein.  Exact  limits 
can  be  obtained  at  the  expense  of  sometime  extensive  computer 
simulation  and  such  simulation  was  used  via  the  computer  programs 
described  in  the  appendices. 

Finally  we  stress  that  the  same  provisos  on  the  models  and 
test  conditions  obtain  as  for  Maximus. 
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2.  CALCULATING  BOUNDS  FOR  A  SINGLE  COMPONENT 

Suppose  that  n  components  of  a  certain  type  are  tested  and 
that  s  prove  to  be  reliable.  If  the  true  reliability 
(probability  of  success)  is  p  then  the  probability  of  the 
above  outcome  is 


PS(1  -  p)n‘S  . 


Let  the  uncertainty  in  p  be  described  by  a  Beta  random  variable  p 
with  probability  density  function 


f(pl  =  p  (1  -  p)  /B(a0,  80),  0  <  p  <  1, 


where 


B(«„  B0)  = 


V1 


P  0  -  p)  dp  . 


This  distribution  is  called  the  prior  distribution.  Using  Bayes' 
theorem  the  posterior  (after  tests)  distribution  is 

f(p)  =  Pa_1  (1  -  P)e_1/B(a,  B)  , 

where  a  =  o0  +  s  and  e  =  80  +  n  -  s. 


Since  the  prior  and  posterior  are  both  Beta,  the  Beta 
distribution  is  said  to  be  conjugate  for  pass/fail  testing.  It  is 
a  flexible  (two  parameter)  distribution  for  expressing  uncertainty 
about  component  reliability  prior  to  testing.  Notice  that  the 
prior  parameter  a0  is  increased  by  the  nunber  of  successes  and 
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3  is  increased  by  the  number  of  failures.  Thus  it  can  be  helpful 
0 

to  relate  a  to  successes  and  3  to  failures  as  in  Appendix  A. 3. 

Let  pa  be  such  that  Prob(p  s  pQ)  =  a, 
i.e.  Prob(p  >  pQ)  =  1  -  a.  Then  is  a  lower  100(1  -  a)% 
Bayesian  limit  for  reliability;  see  Fig.  1. 

Figure  1 


Often  a  uniform  prior  (a0  =  80  =  1 )  is  used  to  express  vague 
prior  knowledge,  whence  the  posterior  probability  density  function 
becomes 

f(p)  =  PS(1  -  p)n'S/B(s  +1,  n  -  s  +  1)  . 

Percentage  points  may  be  obtained  by  interpolation  in  tables  of  the 
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incomplete  beta  function  such  as  Table  16  of  Biometrika,  Vol.  1. 

Exact  limits  can  be  obtained  using  the  computer  program  described 
and  listed  in  Appendix  B. 

Example  1 

Let  s  =  11,  n  =  12  with  a  =  6  =1.  The  posterior 
distribution  of  reliability  is  then  B(12,  2).  Since  the  parameters 
are  integer  Table  16  can  be  used  without  interpolation.  For  a  =  0.1 
the  required  entry  is  that  for  which  v2  =  2  x  12  =  24  and 
v2  =  2x2  =  4.  This  is  0.7322,  which  is  the  lower  90%  limit  for 

reliabi lity. 
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3.  CALCULATING  LOWER  LIMITS  FOR  A  SERIES  SYSTEM  OF  NONREPEATED 
INDEPENDENT  COMPONENTS 

Let  there  be  k  components  in  the  series  system  and  suppose 
that  the  posterior  distributions  for  reliability  are  independent 
Beta  distributions  B(a.j,  g^);  i  =  1,  2,  ....  k. 

The  components  are  numbered  1,  2,  ....  k  so  that 

“i  +  8i  s  “2  +  i  ak  +  sk‘  Tlle  reduction  rules  are 

described  for  k  =  2  and  extend  immediately  for  general  k. 

Rule  (i):  If  >  <x2  +  B2  then  take  a$  =  a2  and 
+  ^1 

a  +  3  =  -  x  (a  +  3  ).  Figure  2  depicts  the  scaling 

s  s  ax  '22' 

down  procedure. 


Figure  2 


a,  +  3  -*■ 

1  i 


a  +  B 
i 


(a2  +  B2)  =  «s  -  Bs 


a 


l 


+  B 


2 


«  -*■  a  =  a 

2  2  s 


Rule  (ii):  If  ax  <  a2  +  B2  take  as  +  &s  =  ai  +  and 


a  a 
1  2 


as  =  a — +~~J~  *  Scaling  down  is  shown  in  Figure  3, 


-2  *2 
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Figure  3 


“x  +  6x 


“l  +  8x  =  “s  +  8S 


+  8 


2 


a 


2 


ct  +3 
2  P2 


a 


S 


When  <*!  =  a2  +  e2  then  a%  +  B$  =  a1  +  B3  and  a$  =  a2. 
In  this  case  the  reduction  is  exact  (see  A. 3). 


Example  2 

Let  k  =  3;  =  28,  Bx  =2;  a2  =  21 ;  B2  =  20,  a3  =  20,  B3  =  2. 

Since  >  a2  +  B2  take  a  +  8  =  (ax  +  81)(az  +  B2 )/aj  =  (30  x  24)/28 

=  25.7143  and  a  =  a2  =21.  Rewrite  a  as  a2  and  8  as  b2  . 

Finally,  since  ~t2  (=21)  <  a3  +  83  (=22),  the  system 
as  =  (21  X  20 )/22  =  19.0909  and  6S  =  6.6234. 

The  lower  90%  limit  for  system  reliability  obtained  by  simulating 
10,000  values  from  the  exact  posterior  distribution  is  0.627.  The 
limit  given  by  the  above  approximating  beta  distribution  is  0.629. 

Exponential  Times  to  Failure 

Suppose  that  a  given  component  has  exponential  failure  rate  X. 

Then,  given  X,  the  time  to  failure  x  has  probability  density 
function 


* 
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f(x)  =  xe  Xx,  x  >  0,  X  >  0  . 

Let  the  uncertainty  in  X  be  described  by  a  gamma  distribution 
with  probability  density  function 

n.-i  “Tox 

g(x)  =  t0(t0x)  e  /r(n0)  » 

t0  >  0,  n0  >  0.  The  likelihood  function  when  n  components  have 
failed  in  a  total  time  on  test  t  is 

*  =  xVxt  . 

The  posterior  distribution  of  the  failure  rate  is  then 
g(X)  =  T(TX}rle‘TV(n)  , 

when  n  =  n0  +  n  and  t  =  xo  +  t.  Again  we  have  a  flexible  prior 
distribution  which  is  conjugate  for  exponential  testing.  Now,  without 
loss  of  generality,  we  may  take  the  unit  of  time  to  be  the  mission 
time  whence  the  reliability  of  the  component  is 


It  follows  that  the  distribution  of  p  has  probability  density  function 
f(p)  =  Tn(-tnp)T1”1pT"1/r(n)  • 

This  is  precisely  the  distribution  of  the  product  of  n  independent 
B(t,  1)  variates.  Thus  such  a  component  in  a  structure  is 
equivalent  to  n  independent  series  components  each  with  a  B(t,  1) 
distribution.  If  n  is  an  integer  than  the  reduction  methods 
given  previously  for  Beta  variates  give  a  single  approximating  Beta 


and  b  =  {(t  +  1)n  -  Tn}/(t  +  1)n_1. 
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distribution  with  a  =  xn/(T  +  I)11"1 
Even  if  n  is  noninteger  the  above  approximation  is  still  very  good. 
Note  that  in  Maximus  noninteger  sample  sizes  and  successes  are  common 
under  reduction. 

Example  3 

Let  o=3,  t  =  40.  Then  a  =  403/41 2  =  38.07  and 
B  =  (41 3  -  403 >/41 2  =  2.93.  The  exact  lower  90%  limit  for 
reliability  is  exp  (-x(.l)}  ,  where  x(.1)  is  the  10%  point 
of  the  distribution  of  X.  Since  it  is  well  known  that  2xt 
has  a  chi-squared  distribution  with  2n  degrees  of  freedom  we  have 
X( .1 )  =  x|(-1  )/2  x  40  =  0.1331  giving  e'^’1  *  =  0.875.  Using 
the  incomplete  beta  function  with  a  =  38.07  and  b  =  2.93  we  get 
the  corresponding  approximate  lower  90%  limit  of  0.864. 


i 
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4.  CALCULATING  LOWER  LIMITS  FOR  A  PARALLEL  SYSTEM  OF  NONREPEATED 
INDEPENDENT  COMPONENTS 

Let  there  be  k  components  in  the  parallel  system  and  suppose 
that  the  posterior  distributions  of  component  reliabilities  are,  as 
before,  independent  B(a.,  8^  distributions;  i  =  1,  2,  ....  k. 

The  reliability  of  the  parallel  system  for  given  component 
reliabilities  px,  p2,  ....  pk  is 

P  =  1  -  nk( 1  -  p.) 
i=i 

=  1  -  n  q  . ,  where  q.  =  1  -  . 

i=i  1  1  1 

Thus  qs  =  1  -  ps  =  nk  qi  . 

If  p^  -  B(a..,  8^  then  qi  -  B(8.j,  ai )  and  the  results  for  series 
systems  can  be  used  here  by  interchange  of  parameters  and  of 
reliabilities  by  failure  probabilities.  Also  one  needs  to  scale 
up  rather  than  scale  down,  as  is  the  case  for  series  sytems.  The 
following  example  illustrates  the  procedure. 


Example  4 

k  =  3; 

al 

=  6, 

9,  e3  =  1. 

Component : 

1 

2 

3 

a  +  8  : 

8 

8 

10 

(a)  : 

(6) 

(6) 

(9) 

8  : 

2 

2 

1 

?■ 

A 

/ 

00 

If 

IS 

10 

(30)  (9) 

2  1 


-  13  - 
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The  rule  for  scaling  up  when  sample  sizes  are  unequal  is  'scale  so  that 
“s  +  ^s  a  nr'n''munl,•  Thus  we  have 

32  10 

(30)  (9) 

2  1 

«s  +  @s  =  16  x  10  =  160 
as  =  (159) 

*S  =  1 

5000  simulated  values  of  system  reliability  gave  a  lower  95% 
limit  of  0.975.  The  lower  95%  limit  given  by  the  approximating 
B(159,  1)  variate  is  0.980. 


i 
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5.  CALCULATING  LOWER  LIMITS  FOR  SERIES  AND  PARALLEL  SYSTEMS  OF 
INDEPENDENT  REPEATED  COMPONENTS 

5(a)  Series  Systems 

Suppose  that  the  posterior  distribution  for  the  reliability  of  a 

component  is  B(a,  8)  and  that  k  such  components  are  connected  serially. 

1/ 

The  reliability  of  the  series  system  is  p  ,  where  p  is  the 
individual  component  reliability.  A  good  approximating  beta  distribution 
has  parameters  as,  8S  given  by 

a(a+1 )...(n+k-1  ) _ 

(a+8)(a+e+1  )...(a+B+k-1 


The  derivation  of  these  expressions  for  ag,  bs  is  given  in  A. 4. 

1/ 

Of  course  it  is  easy  to  obtain  exact  lower  limits  for  p  since 
they  are  simply  the  corresponding  limits  for  p  raised  to  the  power  k. 
However,  when  such  a  structure  is  embedded  in  a  more  complex  arrangement, 
with  components  of  different  types,  we  do  require  a  reduction  to  a 
single  beta  distribution.  The  following  examples  show  that  the  above 

reduction  works  very  well.  It  is  based  on  the  adapted  sequential 
procedure  plus  the  method  of  moments  -  see  A. 4. 

Example  5 

k  =  2;  a  *  30,  8  =  2. 

=  14.762 


/ .  30  x 

V  '  ITT 


Ij}  ‘  1 


2. 
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Thus  B(14.762,  2)  is  the  reduced  beta  variate  for  p2.  The  exact 
lower  90%  limit  for  p2  from  Table  16  of  Biometrika,  Vol.  1  is 
(0.88 023)2  =  0.7748.  The  approximate  limit  given  by  the  above 
reduced  beta  variate  is  0.7747,  in  excellent  agreement  with  the  exact 
result.  Maximus  advocates  equal  splitting  for  repeated  components 

in  series.  The  Lindstrom  and  Madden  method  is  then  applied  as  if  the 
components  were  of  different  types.  The  reduced  beta  variate  is  then 
8(14.06,  1.94),  using  the  analagous  procedures  of  Chapter  3.  The 
approximate  lower  902  limit  for  p2  is  0.770,  which  is  not  as  good  as 
the  adapted  sequential  approach. 


Example  6 


K  =  4;  a  =  30,  8=2 


a 


s 


30  x  31  x  32  x  33V 
32  x  33  x  34  x  35/ 


7.154, 


e 


s 


=  2. 


The  exact  lower  902  limit  is  (0. 88023)"  =  0.6003  and  the  approximate  limit 
is  0.5999,  again  in  excellent  agreement.  Equal  splitting  of  parameters 
gives  0.5770. 


5(b)  Parallel  Systems 

The  formula  for  repeated  components  in  parallel  is  given  in  A.4(ii). 
We  have 


a 


s 


a 


6 


s 


e(e  +  D...(b  +Jc  -  i)  /l'1  . 

(a+@)  (a+B+1 ) . . .  (a+B+k-1 )/ 
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Exact :  The  lower  95%  limit  for  system  reliability  is  1  -  the  upper  95% 

limit  for  failure  probability.  For  q  this  limit  is  0.4295 
and  therefore  exact  lower  95%  limit  for  reliability  is 
1  -  0.0326  =  0.967. 

Approximate:  Using  the  approximating  B(8,  .0563)  variate  the  lower 


95%  limit  is  0.960 
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6.  TWO  COMPLEX  STRUCTURES 


Structure  1 


(1‘,  1"),  (I1,  1"),  (3‘,  3")  signify  replicates  of  components 

1,  2  and  3  respecti vely. 

Beta  Parameters 
Component  a  3 

1  29  1 

2  10  1 

3  9  1 

4  8  2 

5  8  2 

6  8  2 

Using  the  reduction  methods  described  previously 
Ml  =  (1*  *  1")  **?  6(14.5,  t);  (27/2")  =>  B(5,  1);  (37/  3")  B(9, 

M2  ->  B(45.83,  0.167),  M3  B(248,  2). 


0.167) 


Ml  *  M2  *  M3  ->  B(14.5,  1.182). 
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Approximate  lower  95%  limit  =  0.792 
Exact  lower  95%  limit  =  0.809. 

(Based  on  5000  simulated  system  values). 

Structure  2 


/  « 

2  —  a  — 


4 


(2',  2")  are  replicates  of  component  2. 


Beta  Parameters 


Gamma  Parameters 


Component 


40 


25 


Converting  the  Gamma  distributions  to  Series  Beta  distributions  using  the 
results  of  Chapter  3  we  get  for  component  1  a  8(38.025,  1.985)  variate 
and  component  4  reduces  to  a  8(22 .1 1 4,  2.886)  variate.  The  final 
system  approximating  Beta  distribution  has  parameters  38.025  and 
2.151.  This  yields  an  approximate  lower  95%  limit  of  0.8806. 
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5000  simulated  values  of  system  reliability  using  the  computer 
program  gave  a  corresponding  lower  limit  of  0.8859. 
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APPENDIX  A 

An  Alternative  to  Lindstrom  and  Madden 

The  Lindstrom  and  Madden  method  is  recommended  for  series  systems  by 
Maximus.  We  first  describe  an  exact  sequential  procedure  for  series 
systems  and  then  show  how  it  can  be  adapted  to  deal  with  arbitrary  but 
fixed  test  sample  sizes.  The  proof  of  an  inti  resting  property  of  the 
adapted  procedure  is  then  given.  Some  comparisons  with  results  given 
by  Lindstrom  and  Madden  and  also  some  other  approximate  procedures  are 
presented.  Finally  it  is  shown  how  the  modified  sequential  procedure 
can  be  converted  to  the  analagous  Bayesian  procedure  given  in  Ch.  3. 

A.  1  An  Exact  Sequential  Procedure 

Let  there  be  k  independent  components  in  series.  The  component 
types  may  be  arbitrarily  labelled  1,  2,  ....  k.  Suppose  ni 
components  of  type  1  are  tested  and  that  si  prove  to  be  reliable. 

Take  the  sample  size  for  component  2  to  be  n2  =  s1  and  suppose  that  s2 
are  reliable.  Continuing  in  this  way  the  sample  size  for  the  kth 
component  is  nk  =  ski  and  let  sk  of  these  be  reliable.  Then  in 
random  sampling  sk  is  Binomially  distributed  with  nx  trials  and 

If 

parameter  p  =  fl^p...  Write 
sk  -  Bin(n1,  p). 

Proof  that  sk  ~  Bin(n1,  p) 

Let  k  =  2.  Given  n2  =  s2  the  probability  generating  function 

for  s„  is 
2 

S 

G2(z|n2  =  sx)  =  (q2  +  p2z)  ,  q2  =  1  -  p2  . 


* 
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The  unconditional  probability  generating  function  is 


G2(z)  =  l  (q2  +  p2z) 


si=0 


Pi  ^  -  Pi> 


nrsi 


[(q2  +  p2z)p2]  1  (1  -  Pj) 


n,-s, 


=  C(q2  +  P2z)p2  +  (1  -  p2)] 

=  CO  '  PXP2)  +  PiPiZ] 

ni 

=  (q  +  p)  .  p  =  PiP2,  q  =  i  -  p  • 

The  result  for  k  components  in  series  follows  by  induction. 


A. 2  A  Forced  Sequential  Procedure 

Now  let  n2,  n2,  ....  n^  be  arbitrary  but  fixed  sample  sizes  and  let 
s2,  s2,  ....  s^  be  the  corresponding  numbers  of  reliable  components 
obtained  in  tests.  Without  loss  of  generality  label  the  components  so 
that  n2  s  n2  j  . ..  5  n^  .  For  simplicity  let  k  =  2  initially.  Two 
cases  need  to  be  considered, 

(i )  s2  *  n2  and  (ii )  s  <  n?  . 

The  rule  to  'force'  the  sequential  procedure  for  case  (i)  is  -  sample  at 
random  and  without  replacement  from  the  test  results  for  component  1  and 
stop  when  n2  successes  have  been  obtained.  The  average  sample  size 


i 

i 
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to  achieve  n2  successes  is  n2  x  (n  +  1  )/(sx  +  1 ).  This  is  taken  to  be 
the  system  sample  size  and  the  nunber  of  reliable  systems  is  taken  to  be  s2. 
For  case  (ii)  choose  a  random  sample  of  size  sx  without  replacement  from  the 
test  results  for  component  2.  The  average  or  expected  number  of  reliable 
items  is  s2  x  s2/n2  and  this  is  taken  to  be  the  number  of  reliable  systems 
in  a  sample  of  size  n  .  In  both  cases  we  assume  that  the  results  are 
binomial  test  results  and  lower  confidence  limits  for  system  reliability  can 
be  obtained  from  binomial  reliability  tables,  interpolating  as  necessary 
for  noninteger  sample  sizes  and  succcesses  (number  reliable).  For  the 
above  procedure  both  the  sample  size  N  and  the  number  reliable  in  the 
sample  S  are  random  variables  connected  by  the  interesting  result 

E(S)  =  p  E(N), 

where  p  =  p2p2.  The  procedure  extends  easily  to  general  k. 


Proof  that  E(S)  =  p  E(N) 


E(N)  =  n2  l1  p(S  =  sx )  +  nx  p(Sx  <  s2) 


"i  fni+1]K  si,  "i-Sl  *  '  n  Sj  n 

l-  S~rr  S  pi  (1  •  +  n!  I  s  p2  d-Pj) 

i  nz  1  Jl  sx=0  i 
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Using  the  combinatorial  identity 


X 

fx-r 

x-1 ' 

y  j 

!  y  J 

+ 

y-l 

we  can  write 


r 

>x=n2 


fn  +r 

S1+1,  "i-Si 

ni]  Sl,  "l‘Sl  1 

V 

Pi  (1-PX)  =  l 

Sl”n2 

[sj  Pi  ^-Pi5  -  1 

"2 

Pi 


Thus 


E(N)  = 


n,  n1 

—  I 

D  L 

Pl  Sl=n2 


V 

S1  „  1 
Pi  o  -  px)  -  n2! 

V 

s,  , 

n 

\ 

2 

n,  n  -n+1 

‘  /j  \  •*-  *■ 


Pi  ‘  Pj 


n2_1 

"l  I 

V° 


P,  1  (1  -  Pl) 


A  second  application  of  the  combinatorial  identity  to  the  last  term  on  the 
right  hand  side  of  E(N)  gives 


n  nx 
E(N)  =  -i  l 

1  si=n2 


"l 

SJ 


n,  -s. 


"i/,  >  i  ■'i 

Px  0  -  Pl)  -  n 


n2-1 


Px  (1  -  Px) 


ni'n2+1 


("a'1 

n  -1 

+  n  i 
i 

l 

v° 

2 

*1-1 


Px  (1  -  Pl)  ♦ 


V1 

n2-1 


V1.,  ini‘n* 

Pi  O  -  P, ) 


n2  iii 

=  —  y 

p  L 

1  Sl=n2 


S1  „  t  s 

P,  O-p.)  *-p7  sj0  > 


s  n 

Px  1  (1  -  Pl) 


=  E(S)/(Plp2) 


E(S)  =  p  E(N)  . 
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A  Comparative  Example 

k  =  3;  r\1  =  30,  s2  =  28;  n2  =  24,  s2  =  22;  n3  =  20,  s2  =  19. 
Reducing  components  1  and  2  using  case  (i)  gives  n  =  25.714, 
s  =  22.  Combining  these  values  with  those  for  component  3,  again 
using  (i),  gives  the  final  values  n  =  23.377,  s  =  19.  Linear 
interpolation  in  the  tables  of  Cook,  Lee  and  Vanderbeck  (1964)  gives  a 
lower  90%  confidence  limit  of  0.670  for  system  reliability. 

The  Lindstrom  and  Madden  method  gives  n  =  20  and  s  =  16.256 
yielding  a  lower  90%  confidence  limit  of  0.654. 

There  is  no  unique  exact  lower  limit  but  other  methods  suggest 
that  exact  procedures  would  give  values  between  0.68  and  0.69  -  see 
Example  1  of  Winterbottom  (1984). 

Note  that  case  (ii)  is  the  Lindstrom  and  Madden  procedure. 

Overall  the  sequential  approach  is  less  conservative  than  Lindstrom  and 
Madden  and  merits  further  study. 

A. 3  An  Analagous  Bayesian  Sequential  Procedure 

Consider  k  independent  beta  variates  B(a.j,  e^,  i  =  1,  2,  ....  k. 

Suppose  that  i  <*2  +  B2  i  ...  i  ak  +  Bk  and  that  ai  =  a2  +  $2, 

....  a|(_1  =  +  Bk.  Then  the  distribution  of  the  product  of  these 

variables  is  exactly  Beta  B(a,  ,  T  b.-). 

K  i  =  i  1 

Now  (see  Chapter  2),  equating  a  with  s,  a  +  8  with  n 
and  regarding  reliability  as  a  random  variable  with  a,  6  fixed,  the 
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result 

E(S)  =  p  E(N) - >  a  =  E(P)(a  +  B).  This  suggests  the  'forced' 

sequential  procedure  of  Chapter  3  when  the  g^,  Bi  are  arbitrary. 

A, 4  Repeated  Components  in  Series  or  Parallel 
(i)  Series  Systems 

Let  k  =  2  so  that  we  require  an  approximating  distribution  for  p2. 
Suppose  that  the  distribution  of  p  is  B(a,  b). 


Split  a,  @  according  to  the  following  sequential  scheme 
x 

y  =  y 
2 

where  x  +  y  =  a  +  8  and  y  +  z  =  a.  Subtraction  gives  x  -  z  =  e. 
Knowledge  of  y  is  not  required  since  we  take  as  +  Bs  =  x  and  =  z. 
Thus  far  x  and  z  are  not  uniquely  determined.  By  the  method  of 
moments  equate  z/x,  the  mean  of  the  approximating  beta  distribution,  to 
the  mean  of  ps,  i.e. 


z/x 


g(a+l ) 
(g+6)(g+S+l ) 


Using  x  -  z  =  8,  x  =  a$  +  B$.  z  =  a$  and  solving  gives 


a 


s 


q(q+l ) 

(g+8)(g+B+1 ), 


and  8$  =  B. 


The  corresponding  expressions  for  general  k 
q{q+1  ) . . .  (a+k-1 ) _ 


8 


1  - 


(g+S)  (a+B+1 ) . . .  (g+S+K-1 ) 


are 

-l 

-  1 


and  es  =  8. 

(ii)  Parallel  Systems 

Since  failure  probability  q  =  1  -pall  that  is  necessary  to  obtain 
corresponding  results  for  parallel  systems  is  to  interchange  g,  B  in  the 
above  results  for  series  systems. 

Thus 

“s  =  “  r 

amd  B  -  a  -  1  -  -^±1 )-,  (B+k-1) _ 

S  (g+B)(g+6+1 ).. . (g+B+k-1 ) 


APPENDIX  B 


B. 1  Directed  graphs 

It  i s  possible  to  represent  a  control  system  bv  a  directed  graph 
where  the  nodes  are  the  components  In  the  system  and  where  the  directed 
lines  represent  the  flow  of  control  signals. 

To  cater  for  branching  at  the  start  of  a  graph  it  is  convenient  to 
introduce  an  extra  node,  called  the  source,  before  the  first  node. 
Thus  a  source  node  has  no  inputs.  Similarly,  to  deal  with  multiple 
returns  at  the  end  of  the  graph,  an  extra  node  with  no  outputs,  called 
the  sink,  is  added  there.  The  graph  then  has  a  unique  entry  and  a 
unique  exit  and  the  reliability  of  the  whole  graph  is  simply  the 
reliability  of  the  connection  between  the  source  and  the  sink. 

In  order  to  describe  the  graph  the  nodes  are  numbered  sequentially, 
starting  with  zero  for  the  source  and  reserving  the  highest  number  for 
the  sink.  It  Is  then  possible  to  store  the  connection  Information  by 
forming,  for  each  node,  the  set  of  nodes  which  directly  receive  Its 
outputs.  A  convenient  representation  for  this  Is  as  a  linked  list  of 
records  called  KIN  where  each  record  has  Integer  fields  FATHER,  SIZE  and 
a  field  SONS,  consisting  of  an  integer  set,  such  that  the  node  FATHER 
has  SIZE  outputs  to  nodes  called  sons,  the  numerical  values  of  which  are 
held  as  elements  of  the  set  SONS. 

If  there  is  a  linkage  through  a  succession  of  nodes  from  A  to  B 
then  the  corresponding  set  of  node  numbers  Is  called  a  path.  Clearly 
there  may  be  more  than  one  path  between  two  nodes  and  any  node  may  be  a 
member  of  several  paths.  Two  paths  are  said  to  be  independent  If  their 
intersection  is  the  null  set  otherwise  they  are  dependent.  Thus 
dependent  paths  have  nodes  in  common  and  Independent  paths  do  not.  A 
path  which  contains  both  the  source  and  the  sink  Is  cal  led  a  complete 
path,  otherwise  It  Is  called  an  Incomplete  path.  The  type  of  control 
system  which  is  considered  In  this  report  consists  only  of  complete 
paths . 

From  the  Information  in  the  list  KIN,  all  distinct  and  complete 
paths  may  be  found.  These  paths  may  be  manipulated  to  give  the 
algebraic  expression  for  the  overall  reliability  of  the  graph  and 
finally  this  expression  may  be  evaluated  numerically. 

This  process  Is  most  suitable  for  embodiment  In  a  computer  program 
which  Interrogates  the  user  for  the  son  nodes  of  each  father  node, 
requesting  also  reliability  Information  for  each  node  and  which  then 
proceeds  to  perform  the  appropriate  manipulations  on  these  data. 
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B.2  Control  system  reliability 

The  algebra  for  dealing  with  the  reliability  of  combinations  of 
components  in  a  control  system  is  based  upon  the  two  results 

(a)  If  nodes  A  and  B  are  directly  connected  in  series  then  the 
reliability  of  the  combination  Is  r* r»  where  r*  and  r»  are  the 
reliabilities  of  A  and  B, . 

(b)  if  nodes  A  and  B  are  connected  In  parallel  then  the  reliability 

of  the  combination  Is  1  -  (  1  -  r*  )  <  1  -  r»  ) . 

These  results  are  readily  extended  to  more  complex  cases  of 

ser les/para 1 1 e 1  connections. 

If  all  control  systems  were  composed  of  ser les/para 1 1 e 1 
combinations  of  components  then  the  rules  (a>  and  (b)  would  be 

sufficient  to  determine  the  overall  reliability  of  the  system.  Now  the 
graph  In  flg.l  Is  an  example  of  a  control  system  where  It  Is  not 
possible  to  discover  a  pair  of  nodes  which  are  either  In  series  or  In 
parallel.  It  follows  that  some  other  technique  Is  required  to  assess 
the  overall  reliability  of  a  system  in  which  the  nodes  are  connected  in 
an  arbitrary  fashion. 

The  new  technique  for  the  analysis  of  general  control  systems 
starts  with  the  determination  of  all  distinct  paths  from  the  source. 
This  Is  achieved  by  using  a  linked  list  of  records  called  ROUTE,  each 
record  contains  a  field  PATH  consisting  of  the  set  of  nodes  In  a  path 

and  a  field  LASTIN  which  holds  the  last  node  that  had  been  added  to 

PATH.  Initially  ROUTE  consists  of  a  single  record  with  PATH  =  101  to 
represent  the  source  node  and  with  LASTIN  =  0. 

The  Information  in  KIN  is  then  accessed  to  extend  the  list  in  ROUTE 
as  foil ows  : - 

For  each  record  in  ROUTE,  consider  LASTIN  to  be  a  father  node 

having  the  sons  st..s.(..,  then  generate  and  add  to  ROUTE  a  total  of 

(SIZE  -  1)  replicate  records  to  give  SIZE  copies.  If  these  copies  have 
PATH  fields  denoted  by  p,  ..  p. , , .  and  LASTIN  fields  L,  .  L. , , . ,  then  for 
k  *  1  to  SIZE  add  s*  to  pk  and  set  L»  «  s,  .  If  a  son  to  be  added  is 
the  sink,  then  the  assembly  of  that  particular  record  Is  terminated. 
This  process  Is  repeated  until  each  record  in  ROUTE  has  been  assembled. 

All  paths  in  the  systems  under  consideration  are  complete  and  It 
follows  that  any  such  system  will  exhibit  overall  failure  only  If  every 

path  falls.  This  means  that  the  paths  must  be  considered,  In  some 

sense,  to  be  connected  in  parallel  between  the  source  and  the  sink. 
However,  this  concept  needs  to  be  modified  by  taking  account  of  the 
dependences  between  paths. 

To  describe  the  algebra  required  for  path  manipulation  let  P  be  a 
complete  path  and  let  (P)  denote  the  associated  reliability  expression 
obtained  as  the  product  of  the  reliabilities  of  the  nodes  In  P.  Then 
for  a  control  system  without  branches,  there  Is  only  one  path  and  <P> 


is  the  system  reliability.  Also  if  P  and  Q  are  any  two  independent 
paths  where  the  nodes  in  P  are  connected  in  series  with  the  nodes  of  Q 
then  the  reliability  of  the  combination  is  the  reliability  of  the  union 
set  for  P  and  Q,  thus  tP>»{Q)  =  (P  U  Q)  .  It  is  convenient  also  to 

use  the  notation  C  P  -  Q  ]  for  the  set  of  elements  which  are  members  of 
P  but  not  members  of  Q. 

Suppose  now  that  Si  and  S3  are  two  complete  paths  and  it  is 
required  to  derive  their  combined  reliability  expression.  Using  T  to 
denote  the  intersection  set  of  Si  with  S2  ,  it  follows  that  the  nodes  in 
T  are  linked  in  series  with  the  parallel  combination  of  the  nodes  in  the 
sets  [  Si  -  T  J  and  C  S2  -  T  1.  The  reliability  of  this  parallel 

combination  is  then 

US,  -  T  J>  ♦  it  S,  -  T  1)  -  it  S,  -  T  1  U  C  Sj  -  T  ]  >  . 

The  connection  now,  in  series,  of  the  nodes  in  T  gives 
(T)  *  U  S,  -  T  1)  *  (T)*U  Si  -  T  ))  -  { T)  *  U  S,  -  T  3  U  [  S*  -  T  3) 
which  is  equivalent  to 

(T  U  I  S,  -  T  3>  +  <T  U  I  S,  -  T  ])  -  {T  U  U  S,  -  T  ]  U  C  S,  -  T  i  )  ) 
and  this  simplifies  to 

CS,  1  *  (S,)  -  (S,  U  5,1  . 

Now  the  expression  1  -  <  1  -  S,  >*(  1  -  S*  )  can  be  taken  to 
yield  S,  +  S*  -  S,  *S3  ,  so  that  if  the  operator  »  is  replaced  with  the 

union  operator  and  the  operators  +,  -  understood  to  act  only  on  the 

equivalent  reliability  expressions  then  the  required  reliability 
expression  is  obtained. 

The  reliability  for  a  further  complete  path  S3  in  conjunction  with 
Si  and  S2  is  then  given  by  associating  S3  with  the  three  components 
above  giving  successively 

iS,  )  +  iSj  )  -  (S,  U  S3  >  , 

(S,  1  *  (S,(  -  (S,  (J  S,)  , 

and  -  (S,U  Sj)  -  tS3  )  +  (S,  U  S3  U  S3  >  . 

The  sum  of  these  reliabilities  can  be  derived  from  the  expression 
1  *  ( 1  -  Si)«(l  -  -  S3 )  if,  as  before,  the  operator  »  is  replaced 

by  the  union  operator. 

It  follows  that,  with  these  interpretations,  rule  (b>  can  be 
applied  to  complete  paths  to  give  their  combined  reliability  expression. 
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The  procedure  to  be  adopted  for  a  general  control  system  then 
becomes  : - 

(1)  find  all  distinct  complete  paths  , 

(ii)  if  there  are  N  such  paths  denoted  by  S,  ,  Sa ,  ..  ,  S„  then 

form  the  expression 

1  -  1  T  <  1  -  S„  ) 

k  =  1 ,  N 


(ill)  evaluate  the  products  in  this  expression  using  the  given 
interpretations  for  the  operators  , 

(iv)  finally,  translate  the  resulting  expression  in  sets  into  the 
corresponding  reliability  expression. 


For  the  system  of  fig.l  the  construction  of  ROUTE  is  as  follows  :- 

tOl  initial  path,  now  add  the  son  of  node  0, 
CO  1]  add  the  three  sons  of  node  1, 
CO  1  21  CO  1  31  CO  1  41  add  the  son  of  node  2, 
CO  1  2  51  CO  1  31  CO  1  41  add  the  sons  of  node  3, 
CO  1  2  51  CO  1  3  51  CO  1  3  61  CO  1  4]  add  the  son  of  node  4, 
CO  1  2  51  CO  i  3  SI  CO  1  3  61  CO  1  4  61  add  the  son  of  node  5, 
CO  1  2  5  71  CO  1  3  5  71  CO  1  3  63  CO  1  4  61  add  the  son  of  node  6, 
CO  1  2  5  71  CO  1  3  5  73  CO  1  3  6  73  CO  1  4  6  71  add  the  son  of  node  7 


to  give  the  four  complete  paths  represented  by  the  sets  i- 

S,  =  CO  1  2  5  7  81,  Sa  =  CO  1  3  5  7  81, 

Sj  =  CO  1  3  6  7  81  and  S,  =  CO  1  4  6  7  81. 


The  reliability  expression  for  the  system  is  now  found  by  evaluating 

(  1  -  S,  )  (  1  -  S*  >  «  S,  U  Sa  -  S,  -  Sa  +1 

=  CO  12357  81  -CO  1257  81  -CO  1357  83+1, 

then  <  i  -  S,  ) <  1  -  Sa ) <  1  -  S,  )  = 

1-  CO  1257  8]  -CO  1357  81  -CO  1367  83  +  CO  12357  81 
+  CO  1  3  5  6  7  81  . 

and  finally,  1  -  <  1  -  S,  ) (  1-S,)<  1  -Sj><  1  -S«>  = 

CO  1257  81+  CO  1357  83+  CO  1367  81+  CO  1467  81 
-CO  12357  83  -CO  13467  81  -CO  13567  81 
-CO  124567  83+  CO  1234567  83. 
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v 


"phfc  corresponding  reliability  expression, 
and  sink  reliabilities  are  both  unity,  is  then 


assuming  that 
found  to  be  : - 


the  source 


fi  (rjr,  ♦  riTj  +  r,r,  ♦  r«  r*  -  r  a  ra  r5  -  r5  r»  r, 
-  fjrsr,  -  r,r,rj  r,  ♦  firir,rjr,)r, 


B.3  Numerical  assessment  of  reliability 

Uhen  the  nodes  in  a  control  system  or  a  program  have  given 
numerical  values  for  their  reliabilities  then  it  is  a  simple  matter  to 
evaluate  the  overall  reliability  by  direct  substitution  in  the 
reliability  expression. 

However,  if  each  node  has  an  assumed  reliability  distribution  then 
it  should  in  principle  be  possible  to  find  the  overall  reliability 
distribution  as  an  analytic  expression.  The  problem  becomes  one  of 
finding  the  theoretical  distribution  for  a  weighted  sum  of  products  of 
variates  with  known  distributions. 

A  simpler  and  more  practical  alternative  to  this  theoretical 
approach  uses  a  Monte  Carlo  method.  Here  the  reliability  expression  is 
evaluated  many  times  using,  for  each  evaluation,  component  reliabilities 
drawn  at  random  from  their  given  distributions.  These  overall 
reliability  values  are  accumulated  in  the  form  of  a  histogram  which  is 
then  taken  to  be  an  approximation  to  the  required  distribution. 

The  fallowing  tabulations  give  the  results  of  the  Monte  Carlo 
assessment  for  the  reliability  of  Structures  1  and  2,  using  5000 
simulated  system  values. 

Reliability  histogram  for  Structure  i 

reliability  frequency 


range 


0.00 

CO 

0.05 

0 

0.05 

to 

0.  10 

0 

0.  10 

to 

0.  15 

0 

0.  15 

to 

0.20 

0 

0.20 

to 

0.25 

0 

0.25 

to 

0.30 

0 

0.30 

to 

0.  35 

0 

0.35 

to 

0.40 

0 

0.  40 

to 

0.  45 

0 

0. 45 

to 

0.50 

0 

0.  50 

to 

0.  55 

0 

0.55 

to 

0.60 

0 

0.60 

to 

0.65 

1 

0.65 

to 

0.70 

0 

0.70 

to 

0.  75 

9 

0.  75 

to 

0.80 

62 

0.80 

to 

0.85 

260 

0.85 

to 

0.90 

903 

0.90 

to 

0.95 

2070 

0.95 

to 

1.00 

1695 
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Mean  reliability  =  0. 92S3 

Variance  estimate  =  0.001949 


Reliability  histogram  for  Structure  2 


reliability  frequency 


range 


0.00 

to 

0.05 

0 

0.05 

to 

0.  10 

0 

0.  10 

to 

0.  15 

0 

0.  15 

to 

0.  20 

0 

0.20 

to 

0.  25 

0 

0.  25 

to 

0.  30 

0 

0.30 

to 

0.35 

0 

0.  35 

to 

0.  40 

0 

0.40 

to 

0.  45 

0 

0.  45 

to 

0.50 

0 

0.50 

to 

0.  55 

0 

0.55 

to 

0.60 

0 

0.60 

to 

0.65 

0 

0.65 

to 

0.70 

1 

0.70 

to 

0.75 

0 

0.75 

to 

0.80 

3 

0.80 

to 

0.85 

66 

0.65 

to 

0.90 

398 

0.90 

to 

0.95 

1739 

0.95 

to 

1.00 

2793 

Mean  reliability  =  0.9474 

Variance  estimate  =  0.001117 


B.4  Program  listing 

The  following  is  the  Pascal  program  for  the  Monte  Carlo  assessment 
of  control  system  reliability. 

program  contro 1  systems  (Input,  output); 

f  This  program  reads  data  which  represents  a  network  of  components 

with  a  start  component  numbered  zero,  a  finish  component  numbere 
intermediate  components  numbered  1..N-*,  and  one  way  links  between  e 
components  forming  various  routes  from  component  zero  to  component  N. 
A  list  of  all  possible  routes  is  built  up  and  then  an  expression  is 
formed  which  represents  the  reliability  of  the  whole  network  in  terms  of 
the  reliabilities  of  the  individual  components. 

The  program  requests  values  for  the  Beta  or  Gamma  parameters  of 
each  node  and  then  uses  a  Monte  Carlo  technique  to  evaluate  the 
reliability  of  the  network  count  <  =  5000  >  times  using  a  pseudo-random 
generator  to  assign  reliability  values  to  the  nodes,  with  each 
reliability  being  derived  from  the  appropriate  Beta  or  Gamma 


distribution.  A  histogram  Is  assembled  for  the  overall  reliability 
distribution.  1 


const 

max  =  50j  t  maximum  component  number  ) 

count  =  5000;  {  number  of  values  for  the  overall 
reliability  histogram  ) 

type 

range  =  0. . max ; 
collection  =  set  of  range; 

ptr  =  "'node; 
node  =  record 

path  ;  collection; 
coeff,  lastin  :  integer; 
next  :  ptr 
end ; 

point  =  “  f ami  1 y ; 
family  =  record 

father,  size  :  range; 
sons  :  collection; 
kin  :  point 
end ; 

params  =  array  [range, 1.. 23  of  real; 
var 

route,  q  :  ptr; 

result,  mean,  varn  :  real; 

reliability  ;  array  (range!  of  real; 

source,  sink  :range; 

links  :  point; 

betapars  ;  params; 

x,  y,  z,  1,  k,  nofO  : Integer ; 

nodellst  ;  array  (range!  of  -1..1; 

(  a  node  list  element  =  -1  If  the  node  has  not  yet  been  mentioned 

=  0  if  the  node  has  been  given  as  an  Input 

=  1  If  the  node  outputs  have  been  supplied  ) 

hlsto  :  array  (0..191  of  Integer;  (  Array  to  hold  the  results  of  ) 

(  a  Monte  Carlo  simulation  for  ) 

{  reliability  as  a  histogram.  ) 

a  :  array  (1..32!  of  real;  <  The  arrays  a,  t,  h,  d  are  required  ) 

t,h  :  array  (1..31!  of  real;  (  for  the  generation  of  variates  from  ) 

d  ;  array  C6..47!  of  real;  (  the  normal  or  Gamma  distributions.  > 
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function  randomtreal; 

{  This  random  number  generator  is  based  upon  N.P.L.  Report  DITC  6/82  by 
B. A. Ulchmann  and  I. D. Hill.  The  global  variables  x,  y,  z  are  to  be 
given  initial  random  integer  values  which  should  be  less  than  30000.  ) 


var 

w,  ran  :  real; 


begin 

repeat 

x  :=  171 »<  x  mod  177) 
if  x<0  then  x  :=  x  + 
y  :=  172* (y  mod  176) 
if  y<0  then  y  :=  y  + 
z  :=  170* ( z  mod  178) 


-  2* ( x  div  177); 
30269; 

-  35* (y  div  176)  ; 
30307; 

-  63*  <  z  div  178) ; 


if  z<0  then  z  :=  z  +  30323; 
w  :=  x/30269  +  y/30307  +  z/30323; 
ran  :=  w  -  trunc(w) 


unti 1  ran  <>  0; 
random  :=  ran 


end ; 


(  The  following  procedures  FL  GM  GS  GT  GO  are  from  Computing  vol.12, 
pages  223  to  246  (1974)  by  J.H. Ahrens  and  U. Dieter.  They  are 
efficient  routines  for  computing  random  variates  from  a  normal 
distribution  (FL)  and  from  a  gamma  distribution  (GM  GS  GT  GO).  The 
function  betagamma  returns  a  random  reliability  either  as  a  Beta  variate 
or  as  exp(  -  Gamma  variate).  ) 


procedure  f 1 (var  x  :  real); 
labe  1 

i,  4,  6,  13,  14,  17; 

var 

u,  us,  tt,  w,  y,  aa  :  real; 
s,  1  ;  integer; 

begin 

1:  u  :=  random; 

s  ; =  trunc ( 2*u ) ; 
u  :=  trunc ( 32* ( 2*u  -  s)); 

1  : =  trunc(u) ; 
if  i  <>  0  then 
begin 

us  : =  u  -  i ; 
aa  : =  aC i J ; 

4;  if  us  >  till  then 
begin 

w  : =  (us 
goto  17 
end ; 


tt i 3)»ht i  )  ; 


35 


u 

:=  random; 

w 

:=  u«(aCi  +  n  - 

aa )  ; 

tt 

:=  ( w/2  +  aa)*w; 

if 

us  >  tt  then 

goto  17; 

u 

:=  random; 

if 

us  <  u  then 

begin 

us  :=  random; 

goto  4 
end ; 

tt  : =  u ; 
us  :  =  random; 
goto  6 
end 
e  1  se 
begin 

i  :=  6; 

aa  : =  aC32i ; 

u  :=  2*u ; 

wh i 1 e  u  <  1.0  do 

begin 

aa  : =  aa  +  d  C  i  ]  ; 
i  :=  i  +  1; 

if  i  >  47  then  goto  1; 
u  !=  2*u 


end 

! 

u  : 

=  u  -  1 ; 

13: 

w  : 

=  u*d[ i ) ; 

tt 

:=  (w/2  +  aa) *w ; 

14: 

us 

:=  random; 

if 

us  >  tt  then  goto  17 

u  : =  random ; 
if  us  <  u  then 
begin 

u  :=  random; 
goto  13 
end ; 

tt  ; =  u ; 
goto  14 
end ; 

17:  y  :=  aa  +  w; 

If  s  =  0  then  x  :=  y  else  x  :=  -y 
end ; 


procedure  gm<n  :  integer;  var  x  :  real); 


var 

p  :  rea  i  ; 
i  :  Integer; 

begin 

p  :=  1  ; 

for  i  : =  1  to  n  do 
p  :=  p*randora; 
x  : =  -  ln(p> 
end ; 


procedure  gs<-a  :  real;  var  x  :  real); 
label  1; 
var 

p,  b,  u,  xx  :  real ; 


begin 

b  :=  <2.718281828459  +  a ) /2 . 71828 1828459 ; 

1 :  u  : =  random ; 

p  i =  b»u ; 
if  p  >  1  then 
begin 

xx  :=  -  ln((b  -  p)/a); 

if  random  >  exp(<a  -  l)*ln(xx)>  then  goto  1  else  x  :=  xx 
end 
e  1  se 
begin 

xx  :=  exp( ln(p)/a) ; 

if  random  >  exp(-xx)  then  goto  1  else  x  :=  xx 
end 
end ; 


procedure  gt(a  :  real;  var  x  :  real); 
var 

m  :  Integer; 
t,  z,  y  .-real; 

begin 

m  : =  trunc ( a ) ; 
f  :  =  a  -  m ; 

If  m  =  0  then  y  :=  0  else  gm(m,y); 
if  f  =  0  then  z  :=  0  else  gs(f,z); 
:=  y  +  z 


x 

end 


V 
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procedure  go(a  t  real;  var  x  :  real); 
1  abe  1 

2,  7,  6,  8,  10; 


var 

mu,  v,  sig,  sig2,  u,  d,  b,  u,  s,  ss,  xx  .-real; 
begin 


mu  :  = 

a  -  1 ; 

v  ;  = 

sqr  t ( a ) ; 

s  I  g2 

:=  a  +  1 . 632993161855*v 

sig  : 

=  sqrt(sig2) ; 

w  :  = 

s i g2/mu ; 

d  :  = 

2. 44948974278318*  s i g ; 

b  :  = 

mu  +  d ; 

u  :  = 

random ; 

if  u  <=  0.009572265238289  then  goto  8; 
f  1  (  s  )  ; 

xx  :  =  mu  +  s i g*  s ; 

if  (xx  <  0)  or  (xx  >  b)  then  goto  2; 
u  :=  random; 
ss  : =sqr ( s ) /2 ; 
if  s  >-  0  then  goto  6; 

if  u  <=  1  -  ss»((l  -  2*s/v)*w  -1)  then  goto  10  else  goto 
6:  If  u  <=  1  -  ss*(w-l)  then  goto  10; 

7:  It  ln(u)  >  mu*(l  +  ln(xx/mu))  -  xx  ♦  ss  then  goto  2  else 

8:  s  :=  -  ln(l  -  random); 

xx  : =  b» ( 1  +  s/d)  ; 
u  :=  random; 

if  ln(u)  >  mu»(2  +  ln(xx/mu)  -  xx/b ) +3. 7203284924588  -  b 
then  goto  2; 

10 ;  x  : =  xx 
end ; 


7; 

goto  10; 

- 1 n ( s i g«d/b  > 


w 
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function  be tagamma < a,  b  :  real  )  .-real  : 

t  Procedure  to  compute  a  reliability  from  the  random  number  generator. 
A  beta  variate  is  returned  for  the  reliability  if  a,  b  are  positive  with 
the  variate  proportional  to  the  cumulative  distribution  of  x»*(a-l>  « 

(1  -  x)*»(b-l).  If  a,  b  have  been  entered  as  negative  values  then  this 
indicates  that  the  required  reliability  is  to  be  derived  from  a  Gamma 
distribution  and  the  value  returned  is  exp(-y/tbl>,  where  y  is  a  Gamma 
variate  from  the  cumulative  distribution  of  x**(lal  -  1)  *  exp(-x). 

The  value  y/(bl  is  the  corresponding  variate  from  the  cumulative 
distribution  for  the  expression  x*»((a(  -  i)  *  exp(-lbl*x).  > 

var 

x ,  y ,  c  : rea  1  ; 
begin 

c  : =  abs  ( a )  ; 

if  c  <=  3  then  gt(c,x)  else  go(c,x); 

if  a  <  0  then  betagamma  :=  exp(  -x/abs(b)> 

e  1  se 

begin 

If  b  <=  3  then  gt(b,y)  else  goib.y); 
betagamma  :  =  x/(x+y) 
end 
end  1 


procedure  split  (var  route,  hd  :  ptr;  n  :  range); 

<  Split  route  into  two  lists,  one  called  hd  containing  those  records 
for  which  n  is  a  member  of  their  path  field,  the  other  list  called 
route  contains  all  the  other  records.  ) 

var 

P  :  Ptr; 
begin 

if  route  <>  nil  then 
if  n  in  route". path  then 
begin 

P  :=  route; 
route  : =  p" . next ; 
p".next  :=  hd ; 
hd  :=  p; 

split  (route,  hd,  n) 
end 

else  split  ( route". next ,  hd,  n) 
end ; 
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procedure  display  <s  :  collection); 
var 

i  :  range; 
begin 

for  i  :=  0  to  max  do 
if  i  in  s  then  write  (i  :  1,  ’  ’); 

ur i te  1  n 
end ; 


procedure  show  (route  :  ptr); 
begin 

while  route  <>  nil  do 
with  route"  do 
begin 

write  (coeff  :  5,  ’  : 

write  (lastin  :  5,  '  :  ’>; 

display  ( path ) ; 
route  :=  next 
end 
end ; 


procedure  create  (var  route  :  ptr); 

{  Read  in  all  the  data  for  the  network  and  create  the  corresponding 
linked  list  called  route.  1 

var 

i,  j,  m,  n  :  range; 

P,  r  :  ptr; 
kith  :  point; 
s  :  coll ect 1  on ; 
done,  test  :  boolean; 

begin 

(  first  find  all  connections  in  kith  giving  the  father  node,  the  set 
of  sons  of  the  father  and  the  number  size  of  these  sons  ) 

repeat 
n  :  =  0 ; 

while  nodelistCn]  <>  0  do  n  :=  n  +  i; 
node  listen]  : =  1 ; 
nofO  :=  nofO  -  1; 

writeln  (’Number  of  connections  from  node  n  :  1,  '  ?  ’); 

read  In  < i  )  ; 

new(kith)  ; 

kith~. father  :=  n; 

kith*" .size  :  =  i  ; 

kith"",  sons  :  =  C  ]  ; 


1 

i 

I 


If  1  =  O  then  sink  :=  n; 

If  <n  <>  source)  and  (n  <>  sink)  then 
begin 

ur 1 te 1 n ( *  Re  1  lab  1 1 1 ty  parameters  of  node  ',n  :  t,’  ?  ’); 
wr 1 te 1 n< ’ Enter  the  negative  values  for  Gamma  distribution’); 
r ead I n ( be tapar s  C  n , 1 1 , bet a par s  C  n ,  2  3 ) 
end ; 

if  1  >  0  then 
begin 

for  j  :=  1  to  1  do 
begin 

wrlteln  (’Connection  j  :  1,  ’  from  node  ’,  n  ;  1,  ’  ?  ’> 

read  In  (m ) ; 

kith'". sons  :=  kith". sons  +  Cm); 
if  nodellst(m)  =  -1  then 
begin 

node listCml  :=  0; 
nofO  :=  nofO  +  1 
end 
end 
end ; 

kith". kin  :=  links; 
links  :=  kith 
until  nofO  =  0; 

{  now  use  the  connections  In  1 lnk3  to  find  all  pathways  through  the 
graph  from  source  to  sink  J 

new<  route ) ; 
with  route"  do 
begin 

path  ;=  COl; 
coef  f  s  =  1  ; 
lastln  :=  0; 
next  s  =  nil 
end ; 
repeat 

p  :=  route; 
done  :=  true; 
repeat 

If  p  <>  nil  then 

If  (sink  In  p".path)  or  (p". lastln  <  0)  then  p  :=  p*.next; 
test  s  =  p  =  nil; 

If  not  test  then 

test  :=  not  (  sink  In  p".path)  and  (  p". lastln  >  =  0) 
until  test; 

If  p  <>  nil  then 

If  not  (  sink  in  p".path)  and  (p". lastln  >=  0)  then 
begin 

kith  :=  links; 

while  kith". father  <>  p". lastln  do  kith  :  =  kith". kin; 

s  :=  kith". sons; 

for  J  !«  1  to  kith". size  do 
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if  j  <  kith". size  then 

begin 

m  :  =  0; 

while  not  (  m  i n  s  )  do  m  : =  m  + 
s  :  =  s  -  [ml ; 
new  ( r )  ; 
r~  :=  p"; 

r~.path  :=  r". path  ♦  Cm3; 
if  not  (  m  in  p~.path>  then 
begin 

r  *  .  las  tin  :  =  in ; 
done  : =  false 
end 

else  r".  lastin  :=  -m; 
r". next  :=  route; 
route  :=  r 
end 
e  1  se 
begin 

in  ;  =  0 ; 

while  not  (  m  in  s  )  do  m  : =  m  + 
if  not  (  m  in  p~.path  )  then 
begin 

done  :=  false; 
p" . path  :=  p".path  +  Cm3; 
p“ . 1  as t in  s  =  m 
end 

else  p". lastin  :=  -  m 
end 
end 

until  (  p  =  nil  )  and  done 
end ; 
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procedure  comb  (hd  :  ptr;  var  route  s  ptr); 

(  Convert  the  linked  list  of  paths  into  the  linked  list  of  reliability 
terms.  ) 


p,  q,  r  :  ptr; 
s  :  coll ect ion ; 
done  :  boolean; 

begin 

p  ;=  hd; 

q  :=  route; 

whi 1 e  q  <>  nil  do 

with  q "  do 

begin 

new  ( p" . nex t ) ; 
p  : =  p" . next ; 

p".path  :=  path  +  hd".path; 
p". coeff  ;=  -  coeff  «  hd". coeff; 
p" .  1  as  t in  : =  sink; 
q  :=  next 
end ; 

pA. next  ! =  nil; 
while  hd  <>  nil  do 
begin 

s  : =  hd" . path ; 
q  ;=  route; 
done  ;=  false; 
repeat 

if  q" . path  =  s  then 
begin 

q".  coeff  ;=  q". coeff  +  hd". coeff; 

if  q". coeff  =  0  then 

begin 

if  q  »  route  then  route  :=  q".next 
else  r".next  :=  q".next; 
dispose  (q) 
end ; 

done  ;=  true 
end 
el  se 
begin 
r  :=  q; 
q  :=  q".next 
end 

until  done  or  (q  9  nil); 
if  done  then 
begin 
r : =hd ; 

hd  :=  hd".next; 
dispose  (r> 
end 
e  1  se 


l 


? 

t 
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•  1  begin 

r * . nex t  : =  hd ; 
hd  :  =  hd* .  nex  t ; 
r*.next*.next  :=  nil 
end 
end 
end ; 


procedure  combine  (route  :  ptr); 
var 

p,  hd  :  ptr; 
begin 

hd  :=  route*. next ; 
route*. next  :=  nil; 
while  hd  <>  nil  do 
begin 

p  : =  hd* . next ; 
comb  (hd,  route); 
hd  ;  =  p 
end 
end ; 
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procedure  evaluate  (route  :  ptr;  var  result  :  real); 

(  Compute  the  numerical  expression  for  the  overall  reliability 
var 

p  :  ptr; 

x  :  real ; 

s  :  coll ect ion ; 

1  s  integer; 

begin 

resul t  : =  0. 0; 

p  :=  route; 

whl 1 e  p  <>  nil  do 

begin 

s  :  =  p~ .  path ; 
x  : =  p * . coef  f  ; 

1  :=  0; 

while  s  <>  Cl  do 
begin 

if  l  in  s  then 
begin 

x  :=  x  *  re  1 iabi 1 i tyt i 3  ; 
s  ;=  s  -  (11 
end ; 

I  :=  i  +  1; 

end; 

result  :=  result  +  x; 
p  :=  p~.next 
end 
end ; 


begin 

nodel isttOl  s  =  0; 

for  i  :=  1  to  max  do  nodelisttil  :=  -1; 
nofO  :=  1; 
at  1 J  :=  0; 

at  21  :=  0.03917608550309;  aC3]  ;=  0.07841241273311; 
a(41  .•=  0.11776987457909;  at5J  :=  0.15731068461017; 
a(6]  :=  0.19709908429430;  at7]  :=  0.23720210932878; 
at81  :=  0.27769043982157;  aC9J  :=  0.31863936396437; 
at  101  s*  0.36012989178957;  atlll  ;=  0.40225006532172; 
at  121  s=  0.44509652498551;  aC131  :=  0.48877641111466; 
at  141  s*  0.53340970624127;  aC151  :=  0.57913216225555; 
at  161  :=  0.62609901234641;  atl71  :=  0.67448975019607; 
atl81  s=  0.72451438349236;  atl91  :=  0.77642176114792; 
at 201  :=  0.83051087820539;  aC211  :=  0.88714655901887; 
at 221  :=  0.94678175630104;  aC231  :=  1.00999016924958; 
at 241  :=  1.07751556704027;  at251  ;=  1.15034938037600; 
at 261  :»  1.22985875921658;  at271  ; =  1.31801089730353; 
a (281  :=  1.41779713799625;  a(291  ;=  1.53412054435253; 
a (301  s*  1.67593972277344;  aC311  ;«  1.86273186742164; 
a(321  :=  2.15387469406144; 


at  n  ) 


for  1  : =  1  to  31  do 
begin 

till  :=  < (at  1  +  11  -  at  1  3  )  /2  ♦  a  t  1  3  >  *  ( a  t  1  +  1  3  - 
hill  :=  (ati+13  -  at i 3 ) / (  1  -  till) 
end ; 

d  t6 1  :=  0.26368432217502;  dt73  :=  0.24250845238097; 

dt81  :=  0.22556744380930;  dt93  ;=  0.21163416577204; 


dt  101 

=  0.19992426749317;  dtlll 

=  0.18991075842246 

dt  123 

=  0.18122518100691;  dtl33 

=  0.17360140038056 

dt  143 

=  0.16684190866667;  dtlS3 

=  0.16079672918053 

dt  163 

=  0.15534971747692;  dtl7J 

=  0.15040938382813 

dt  183 

=  0.14590257684509;  dtl93 

=  0.14177003276856 

dt203 

=  0.13796317369537;  dC213 

=  0. 13444176150074 

d  1 22  3 

=  0.13117215026483;  dt233 

=  0.12812596512583 

d  1 24  3 

=  0.12527909006226;  «£253 

=  0.12261088288608 

d  1 26  3 

=  0. 120103SS9656S1 ;  dt273 

=  0. 11774170701949 

d  1 28] 

=  0.11551189226063;  dt291 

=  0.11340234879117 

d  1 30  3 

=  0.11140272044119;  dt313 

=  0.10950385201710 

d  1 32  3 

=  0. 10769761G56476;  dt333 

=  0. 10597677198479 

dt343 

=  0.10433484129317;  dt353 

=  0. 10276601206127 

d  1 36  3 

=  0.10126505151402;  dt373 

=  0.09982723448906 

d  1 38  3 

=  0.09844828202068;  dt393 

=  0.09712430874765 

d  1 40  3 

=  0.09585177768776;  dt413 

=  0.09462746119186 

d  1 42  3 

=  0.09344840710526;  dt433 

=  0.09231190933664 

d  1 44  3 

=  0.09121548217294;  dt453 

=  0.09015683778986 

dt463 

=  0.08913386650005;  dt473 

=  0.08814461935364 

x  ;=  31415; 
y  :=  27188; 

wr 1 te 1 n( ' 1 nput  random  seed,  an  integer  between  0  and  30000’) 
read  1 n ( z) ; 
source  ;=  O; 
sink  :=  max; 

1  inks  : =  nil; 
create  ( route ) ; 
q  ;=  route; 
comb lne ( route )  ; 

show (route);  write  In;  write In ; 

rel iabl 1 ItyCOl  ;=  1.0; 

re  1 labl 1 1 ty t s ink ]  ;=  1.0; 

for  1  : =  0  to  19  do  hlstoti]  ;=  0; 

mean  :=  0.0; 

varn  :=  0.0; 

for  1  :=  1  to  count  do 

begin 

for  k  : =  1  to  slnk-1  do 

re  1 labl 1 1 ty t k 3  :  =  betagamma ( be tapars t k , 1 3 , betapars t k, 2 3 )  ; 
evaluate  (route,  result); 
mean  :=  mean  ♦  result; 
varn  :=  varn  ♦  sqr(result); 
k  :=  trunc(20*resu 1 t ) ; 
hlstotk3  :=  hlstotk)  +  1 
end ; 

wr 1 te ln( ’ 


range 


f  requency’ ) 


-  46  - 


for  1  :  =  0  to  19  do 

wrltelni i»0. 05:7:2, ’  to  ' ,  ( i *  1 ) *0. 05 : 4 : 2, h 1 s to C 13  )  ; 
mean  :=  mean/count; 

varn  :=  (varn  -  count  »  sqr (mean ))/( count  -  1); 
writeln;  writelnt*  Mean  reliability  =’,mean); 
wrlteln(’  Variance  estimate  =’,varn) 
end. 


B.5  Procedures  to  compute  beta  and  gamma  functions 

The  following  Pascal  procedures  are  appended  in  case  numerical 
values  are  required  for  Beta  or  Gamma  functions.  They  do  not  form  part 
of  the  program  for  analysing  the  control  structures. 

function  logama<  x  :  real ) :real ; 

{  compute  log  gamma ( x )  function  > 
var 

xc,  f ,  z  : rea 1 ; 

begin 

xc  : =  x ; 

if  xc  >=  7.0  then  f  :=  0.0 

e  1  se 

begin 

f  !=  1.0; 
z  ;  =  xc ; 

whi 1 e  z  <  7.0  do 
begin 

xc  : =  z  j 
f  :=  f*z; 
z  :=  z  +  1.0 
end ; 

xc  : =  xc  +  1.0; 
f  -  l n ( f ) 
end ; 

z  : =  1 . O/sqr ( xc) ; 

logama  :=  f  +  ( xc-0. 5 ) » 1 n < xc)  -  xc  +  0.918938533  + 

<  C (4.0*z-3.0«z) *z- 14.0) »z  +  420. 0)/< 5040. 0*xc> 

end ; 

function  logbetai  p,  q : rea 1 ) ; rea 1 ; 

{  compute  complete  log  beta  function  ) 
begin 

logbeta  s*  logama(p)  ♦  logama(q)  - 
end ; 


logama ( p+q > 


function  betaln(  x,  p,  q,  zeta 


rea I ;  var  1 f au 1 t 


Integer) : rea I 


1  compute  Incomplete  beta  function  ratio  for  argument  x  between 
0  and  1  and  for  p,  q  positive  ) 

1  abe  1 

3,  4,  5; 

const 

accuracy  =  1.0e-8; 
var 

index  , ok  :  boolean; 

psq,  cx,  pp,  qq,  xx,  term,  ai,  rx,  temp,  beta  :  real; 
ns  :  integer; 

begin 

beta  :=  x; 
ifault  :=  0; 

if  (p  <=  0.0)  or  <q  <=  0.0)  then  ifault  :=  i; 
if  (x  <  0.0)  or  <x  >  1.0)  then  ifault  :=  2; 
ok  : =  ifault  =  0 ; 

if  ok  and  (x  <>  0.0)  and  (x  <>  1.0)  then 
begin 

psq  i=  p  ♦  q; 
cx  :=  1.0  -  x; 
index  :=  p  <  psq*x; 
if  index  then 
begin 

xx  : =  cx ; 
cx  : =  x ; 

PP  s=  q? 
qq  s=  p 
end 
e  1  se 
begin 


XX 

;  =  x; 

pp 

:=  p; 

qq 

q 

end ; 

term 

:=  1.0; 

ai  ;  = 

1.0; 

beta 

:=  1.0; 

ns  ;  = 

trunc(qq  +  cx«psq); 

r  x  s  = 

xx/cx ; 

temp 

:=qq  -  ai; 

if  ns 

=  0  then  rx  :=  xx; 

term 

:=  term»temp*rx/(pp+ai ) 

beta 

;=  beta  +  term; 

temp 

: =  abs ( term  > ; 

if  (temp  <=  accuracy)  and  (temp  <=  accuracy»be ta )  then  goto  5 
ai  :=  al  +  1.0; 
ns  : =  ns  -  1 ; 
if  ns  >=  0  then  goto  3; 
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temp  :=psq; 

psq  : =  psq  +  1 . 0  j 

goto  4; 

5:  beta  :=  be ta*exp ( pp* 1 n< xx )  +  ( qq- 1 . 0 ) *  1 n < cx )  -  2eta)/pp; 

If  Index  then  beta  :=  1.0  -  beta 
end ; 

betain  :=  beta 
end ; 

function  betafn(  x,  p,  q  : real ) : rea i ; 

{  compute  Incomplete  beta  function  ratio,  see  change  below  ) 

{  for  how  to  compute  the  Incomplete  beta  function  Itself.  } 

var 

a,  r  :  rea 1  ; 
j  :  Integer; 

begin 

r  ; =  logbetalp, q) ; 
a  :=  betalnl x , p, q, r , J )  ; 

If  J  =  O  then  betafn  :=  a  1  for  beta  function  replace  :=  a  ) 

1  w i th  : =  a*exp ( r )  ) 

e  1  se 
begin 

wr 1 te 1 n (’ Error  In  Incomplete  beta  function’); 
writelnl’ x=’ , x:4, ’  should  be  between  0  and  1’); 
wr 1 te 1 n ( ’ p= ’ , p: 4, ’  should  be  >  O’); 
wr 1 te 1 n ( ’ q= ’ , q : 4, '  should  be  >  O’) 
end 
end ; 
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